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Abstract 

Background: Existing incidence estimates of heroin use are usually based on one information source. This study 
aims to incorporate more sources to estimate heroin use incidence trends in Spain between 1971 and 2005. 

Methods: A multi-state model was constructed, whereby the initial state "heroin consumer" is followed by 
transition to either "admitted to first treatment" or to "left heroin use" (i.e. permanent cessation or death). Heroin 
use incidence and probabilities of entering first treatment ever were estimated following a back-calculation 
approach. 

Results: The highest heroin use incidence rates in Spain, around 1.5 per 1,000 inhabitants aged 10-44, occurred 
between 1985 and 1990; subdividing by route of administration reveals higher incidences of injection between 
1980 and 1985 (a mean of 0.62 per 1.000) and a peak for non-injectors in 1990 (0.867 per 1,000). 

Conclusions: A simple conceptual model for heroin users' trajectories related to treatment admission, provided a 
broader view of the historical trend of heroin use incidence in Spain. 
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Background 

Recently, there has been increasing interest in ascertain- 
ing illegal drug use incidence for planning and evaluat- 
ing prevention strategies [1]. In the case of heroin use, 
as survey data is not effective [2,3], incidence has been 
estimated from users eventually showing up mostly in 
treatment registers [1,4-9]. This incidence can be re- 
ferred to as "problematic use incidence" and its trend 
could provide a satisfactory overview of problematic her- 
oin use, assuming a constant proportion over total 
incidence. 

The present study is an attempt to take the definition 
of problematic use incidence further by incorporating a 
proportion of individuals that used heroin who never 
show up in the main source under study (usually treat- 
ment). This proportion is based on other kinds of 
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information, aggregated sources, estimates, or assump- 
tions, such as mortality and cessation rates. 

The idea is to study the unobserved entry (or immigra- 
tion) of people to the state of "consumer", based on a 
later first entry to treatment. Since heroin users may exit 
the state of consuming heroin before entering treatment, 
whether due to death or permanent cessation of their 
consumption, a "left heroin use" state is added to the 
model. The situation is thus represented by a set of mu- 
tually related states, in a so-called multi-state model 
with immigration. 

This approach has parallels in studies in the HIV field, 
where a multi-state model was presented describing pro- 
gression of HIV disease from infection to AIDS in sev- 
eral stages [10,11] or Rossi's dynamic "mover-stayer" 
model as a theoretical approach applied to simulate a 
complete drug user "career" [12,13]. 

We shall use a back-calculation type approach to esti- 
mating the incidence. This is similar in spirit to de 
Angelis et al. [14]. However, following the presentation 
of Aalen et al. [10], we have found it useful to display 
the approach as a simple multi-state model. 
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To our knowledge, previous approaches had to assume 
that treatment availability was stable over time [1,4-9]. As 
treatment data is usually only available for a limited period 
of time, a lag time distribution between heroin use onset 
and first treatment ever has been employed to avoid 
underestimating incidences. However, this assumption 
may be too strong if important changes in treatment avail- 
ability had occurred. Therefore, the back-calculation ap- 
proach used in this study avoids this assumption. 

In Spain, heroin use generated important health pro- 
blems in the eighties and early nineties [15], thereafter 
all indicators (mortality, treatment admissions, hospital 
emergencies, surveys, etc.) showed a decreasing trend 
until 2006 [16]. Efforts to calculate incidence of heroin 
use have been and still are considerable to understand 
the overall trend and assess its consequences [1]. To 
date only one previous study has estimated heroin use 
incidence in Spain. However, due to important changes 
of treatment availability, the reliability of the estimation 
was questionable [6]. 

The multi-state model approach may be successfully 
used to estimate heroin use incidence with the available 
Spanish heroin users' treatment data. We use the same 
treatment database as the previous study to assess differ- 
ences in the incidence estimates depending on the ap- 
proach used. 

The objective of the present study was to estimate her- 
oin use incidence in Spain through a multi-state model 
with immigration and assess differences with previous 
study's estimates. 

Methods 

Incidence estimation: the multi-state model with 
immigration 

We will describe a heroin user's trajectory simply using 
the time of transition from state 1 (first heroin use) to 
one of two possible subsequent states: first treatment 
ever (state 2), or leaving heroin use before any treat- 
ment, either by permanent cessation or death (state 3) 
(Figure 1). 

We shall now formulate our back-calculation type 
model. We let t be the calendar time and assume that the 
number of new heroin users per unit of time is Poisson 



distributed. The Poisson assumption is standard and well 
justified on statistical grounds [14]. Moreover, we let h t de- 
note the expected number of people entering state 1 (her- 
oin use) at t. The number p t denotes the probability that a 
given heroin user initiates their first treatment ever at time 
t, given that they were in state 1 at the previous time. We 
assume this probability independent of the era when their 
heroin use began. The number q t denotes the probability 
of an individual leaving heroin use at time f, given that 
they were in state 1 at the previous time. The cause could 
be death or other permanent cessation of heroin use. We 
want to estimate the parameters h t , the expected number 
of new heroin users at time t. The probability q t is 
assumed known. 

We now construct the likelihood function. This is a 
standard Poisson type likelihood, following the approaches 
in the basic back-calculation papers [14]. However, since 
our model is here adapted to the particular data we have, 
we briefly present the necessary formulas. Let AT,y be the 
observed number of individuals that start heroin use in 
year i and enter first treatment in year /'. Let be the 
expected value of Ny A short computation shows that: 



k=i 



Pk ~ W-Pj 



(1) 



We see that y,» is the product of the expected number 
of new heroin users in year i {hj), the probability that a 
given heroin user remains in state 1 from time i to time 
j-1 and the probability of a transition from state 1 to 
state 2 at time It also follows that each Nu is Poisson 
distributed. This gives us the following simple expression 
for the likelihood: 

L = n^-exp(-^) (2) 

ij 

Maximizing this likelihood yields estimates for h t and p t . 
Treatment data 

The Spanish Drug Observatory maintains a drug infor- 
mation system. Its indicator "treatment" is based on data 
from all treatment starts in public and publicly funded 



h t 
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—2^ 
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heroin use 




first treatment ever 
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left heroin use 



Figure 1 Multi-state model diagram. Parameter h, denotes immigration of individuals starting heroin use, p t transition rate entering treatment 
for the first time, and q t transition rate leaving heroin use without entering treatment, at time f. 
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centres. The health coverage is universal and all kinds of 
treatment are considered. In the present study, we 
included 169,257 persons who entered first treatment 
for heroin use from 1991 to 2006 when aged 15 to 54 
(mean 28 years), and who started their heroin use be- 
tween 1971 and 2006 when aged 10 to 44 (mean 21 
years). The database was split into two subsets: people 
who declared injection as the most frequent route for 
heroin use in the last 30 days before first treatment ad- 
mission (29%) and people who declared using routes 
other than injection (68%). Route of administration was 
missing in 3% of persons. 

Assumptions about parameters 

As in all applications of back-calculation where a 
detailed history of individuals is not observed [14], we 
have to make some simplifying assumptions, presented 
below. 

First treatment data for the period 1971-1990 

In the estimation of p t , first treatment data was avail- 
able for t between 1991 and 2006, thus restricting esti- 
mates to this period. For the preceding part of the 
study period (t in 1971-1990), we made an educated 
guess of p t based on general heroin use information in 
Spain. Based on the first appearances of admissions for 
heroin use in the emergency units in Spain in 1982 
[17-19]; we assumed probabilities of entering treatment 
(p t ) as low as 0.01 between 1971 and 1981, as there 
were still no specific treatments available. Thereafter 
we assumed a linear increase to the value estimated for 
the parameter p t in 1991. 

Mortality for heroin users 

Mortality rates for heroin users were only available from 
two local cohort studies covering the period 1985 to 
1999 [20,21] in an area where injecting was the predom- 
inant route of administration [22]. As we did not have 
better approximation, yearly rates from these studies 
(minimum 1.4% in 1985, maximum 6.6% in 1995) were 
extrapolated to the whole country for the corresponding 
year. For the period 1971 to 1984 a smooth increasing 
trend from a mortality rate of 1% to 1.4% was applied. 
Mortality rates from 2000 to 2006 decreased from a rate 
of 1.5% in 1999, to 1% in 2006. In the analysis by route 
of administration the same mortality rate was used for 
injectors, but a constant mortality rate of 1% for non- 
injectors since they have lower risk [23]. 

Degenhardt et al. reported a pooled crude mortality 
rate of 2.09 per 100 person-years and that mortality risk 
was increased among out-of-treatment heroin users [24]. 
In a sensitivity analysis (see analysis section) we consid- 
ered alternative mortality rates obtained by adding 0.01 
to the yearly mortality rates in order to ensure a 



minimum rate of at least 2%. Note that, in the multi- 
state model we are imputing mortality rates before first 
treatment. 

Cessation rates 

Owing to the impossibility of obtaining permanent ces- 
sation rates, we looked for lasting cessation rates from 
long-term cohort studies. As such studies are not avail- 
able in Spain, we considered yearly cessation rates from 
a thorough review which reported a range of 0.02-0.04 
[25]. Our analyses considered these two extreme values. 

Analysis 

For all heroin users and for injectors and non-injectors 
separately, we applied the aforementioned multi-state 
model with the Spanish treatment data and the 
assumed leaving rates to estimate the heroin use inci- 
dence (h t , t ranging from 1971 to 2006). As explained, 
the model also estimates the probability of entering 
first treatment {p t , t in range 1991 to 2006). We consid- 
ered the yearly cessation rate of 0.04 and the non- 
modified mortality rate derived from the local cohort 
studies. Note that the probability of leaving drug use 
without having ever been registered for first treatment 
(q t ) is the sum of the cessation and mortality rates for 
each year from 1971 to 2006. 

In equation 1 when i=j, it means that users began 
treatment in the same year as they started heroin use. 
This gives on the average about half a year of observa- 
tion, and so we must weight fiij by 0.5. 

To assess the fit of the expected incidence values fig 
with their observed values Ny, we have drawn their 
curves stratified by year of heroin use onset (i). 

As results can be dependent on assumptions, a sensi- 
tivity analysis was performed to evaluate the two chosen 
mortality and cessation rates obtaining four combina- 
tions of q t , that are reflected in four curves of estimated 
incidence rates. These combinations were: firstly and as 
a matter of choice, the available mortality rates and a 
yearly cessation rate of 0.04; secondly, the same mortal- 
ity rates and a yearly cessation rate of 0.02; thirdly, the 
same mortality rates modified by adding 0.01 to the rate 
for each year and a yearly cessation rate of 0.04; and fi- 
nally, the modified mortality rates and a yearly cessation 
rate of 0.02. 

Statistical uncertainty was estimated using a bootstrap 
technique with 500 re-samples, where each re-sample 
was made up of two parts: 1) the treatment database 
was re-sampled with replacement and, 2) both our "best 
guesses" for p t in the period from 1971 to 1990, and the 
cessation rate for all years were sampled from gamma 
distributions. The shape and scale parameters were 
derived from the mean and standard deviation, taking 



Sanchez-Niubo et al. BMC Medical Research Methodology 201 3, 1 3:4 
http://www.biomedcentral.eom/1471-2288/13/4 



Page 4 of 9 



the mean as the "best-guess" value, and the standard de- 
viation was established as 0.01. 

The expected number of new heroin users per year 
was obtained, and converted into rates per 1,000 inhabi- 
tants, based on Spanish population yearly census data 
for people aged 10-44 [26]. 

Incidence estimates from the previous study were 
retrieved to compare with the present estimates. Both 
the period of years covered and the census figures were 
the same. 

The software used for the statistical analysis was R ver- 
sion 2.13.0 [27]. The study was approved by the Ethics 
Committee of the Institution with number 2004/1828/1. 

Results 

Estimates of incidence and probability of entering 
treatment 

Applying the multi-state model with immigration yielded 
estimated probabilities of entering treatment for the first 
time (p t ) which exhibited an overall increasing trend 
from 0.08 (95% CI 0.07-0.09) in 1991 to 0.29 (95% CI 
0.23-0.49) in 2005. Incidence estimates of general heroin 
use and by route of administration with 95% confidence 



intervals are plotted in Figure 2. For general heroin use, 
the highest incidences were between 1985 and 1990 with 
rates around 1.5 new heroin users per 1,000 inhabitants 
aged 10-44, followed by a steep decline from 1991 to 
1997, then a more gradual decrease from 1998 (0.24 per 
1,000) to 2005 (0.05 per 1,000). 

In the analysis by route of administration the probabil- 
ity of entering first treatment for heroin users declaring 
injection was higher than for non-injectors between 
1991 and 2000 with a difference of around 0.03. How- 
ever, after 2001 this difference increased progressively 
(Figure 3). Incidence rates for injectors were higher than 
for non-injectors until 1985 and lower thereafter (lower 
graph in Figure 2). For injectors the highest values were 
observed between 1980 and 1985 (a mean of 0.62 per 
1,000) whereas for non-injectors the peak was in 1990 
(0.86 per 1,000). 

Comparing the curves of the expected incidence values 
Hij with their observed values Njj, by year of heroin use 
onset (z), we assessed that the fit was good (Figure 4). 

We could also check this good fit modifying the observed 
values as A/y.j and the expected values as ^y_* , where j-i 
represents the lag time between drug use onset i from 1991 
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t Non-injectors include mostly smoking and snorting 

Figure 2 Estimated heroin use incidence rates*: global (upper graph) and by route of administration** (lower graph). * rates per 1,000 
inhabitants aged 10-44, in Spain, with 95% confidence interval, ** the most frequent route of administration in the last 30 days before entering treatment. 
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Figure 3 Estimated probabilities of entering first treatment ever (p t ) in Spain by route of administration**, with 95% confidence 
interval. ** The most frequent route of administration in the last 30 days before entering treatment. 



to 2004 and first treatment ever /, conditional on treatment 1980 and, although they were lower than present ones in 
starting before 2006 (Figure 5). the 90s, they overtook them in the last few years. 



Sensitivity analysis 

Each combination of mortality and cessation rates pro- 
duced large variations in the eighties when incidence was 
highest (Figure 6). Note that lower cessation and mortality 
rates yielded lower incidence rates, and vice versa. 

Regarding probabilities of entering treatment for the 
first time, from 1991 to 2000 estimates only varied slightiy, 
with a maximum difference of only 0.02 (data not shown). 
From 2001 differences increased progressively reaching a 
final value of 0.2 in 2005 between the lowest and highest 
combinations of cessation and mortality rates (probability 
of 0.22 versus 0.42). 

Confidence intervals for the various estimates overlap, 
except for the combination yielding the lowest estimates 
(results not shown). 

Comparison of incidence estimates 

Incidence estimates from a previous study [6] are shown 
in Figure 6. These estimates had an earlier peak around 



Discussion 

We have established a conceptually simple multi-state 
model to obtain estimates for the incidence rates of heroin 
use and applied it over a long period in Spain. The highest 
incidences were observed from 1980 to 1985 correspond- 
ing to injectors and a peak in 1990 to non-injectors. 

In comparison with previous studies, our estimates are 
wider in scope since by including mortality and other 
permanent cessation into the multi-state model it is pos- 
sible to account for almost all problematic heroin users 
after drug use onset. 

The conceptual model employed in this study focuses 
on the first phase of a heroin user's "career": from heroin 
use initiation to treatment. Other more complete models 
based on a theory of compartmental epidemic models 
over drug user "career" have also been described [12]. 
Adding more states into the model using the data avail- 
able would, however, make estimating incidence too 
complex because heroin use cessation and relapse are 
frequent and difficult to follow up. Knowing first entry 
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Figure 4 Observed and expected number of individuals by heroin onset cohort from 1971 to 2005, and entering first treatment ever 
from 1991 to 2005, in Spain. 



to treatment and first use, we only need to account for 
quitting heroin use before first treatment through 
complete recovery (cessation) or death. 

As a Spanish anthropologist described, in Spain heroin 
use had its first phase in the years 1977 and 1978, when 
the first users became visible, being endemic in the second 
phase between 1979 and 1982, and reaching its zenith be- 
tween 1983 and 1986 leading to the institutionalization of 
the problem [18]. Therefore, trends in overall heroin use 
incidence obtained in the present study seem reasonable 
and consistent with previous knowledge about the Spanish 
heroin epidemic and the HIV- AIDS epidemic [15]. Specif- 
ically, the inflection point observed in 1985, when the inci- 
dence rate of injected heroin use fell below the rate for 
non-injection, is consistent with the trend of decreasing 
HIV incidence among injectors in Spain. However, we ob- 
serve that estimates from the previous study had higher 
incidence figures earlier than the present ones (Figure 6). 
They reflect the fact that the availability of treatment was 
assumed stable throughout the entire period, leading to 
high estimates too soon. 



We found a decreasing trend in the incidence estimates 
for the last years observed, which is probably related to 
the decreasing trend observed in all indicators towards 
the end of the period studied, as mentioned in the intro- 
duction. However, estimates for these last years from the 
previous study became stable overtaking the estimates 
from the present one (Figure 6). This is due to the two 
studies employing different approaches. Equation 1 in the 
present study was formulated assuming that p t , the prob- 
ability of entering first treatment, was independent of the 
era when a person's first heroin use began. Actually, this 
would be not entirely true if lag time between the drug 
use onset and first treatment followed a determined pat- 
tern, as previous studies assumed [1,6]. However, if we ob- 
serve Figure 5, the lag time distribution for the observed 
values 2Vy_i and for the expected values fiy-i (j-i represents 
the lag time), for each year of heroin use onset from 1991 
to 2004 all fitted well. So, to modify the equation 1 includ- 
ing the probability of entering first treatment conditional 
on the initiation of heroin use would be too complex and 
may not have great practical importance. Therefore, bias 
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Figure 5 Observed and expected distributions of lag time between heroin use onset and first treatment ever conditional on treatment 
starting before 2006. Distributions are given by heroin onset cohort from 1991 to 2004, in Spain. 



can be inherent in both the independence assumption and 
assuming a determined pattern of lag time, although the 
direction of such bias cannot be determined. 

As observed by de Angelis [14], results are dependent 
on assumptions. However, the sensitivity analysis showed 
that the different incidence curves generated by varying 
the cessation and mortality rates had similar shapes 
(i.e. trends) although different levels, suggesting that the 
model's estimates were stable. Moreover, we observed 
that the confidence intervals of incidence figures esti- 
mated using non-modified mortality rates and a cessa- 
tion rate of 0.04 completely contain those of the other 
three estimates. Thus the chosen incidence estimates do 
not differ significantly from the other three estimates, 
which resulted from varying the rates involved. 

Concerning the assumptions made about the model 
parameters, such as 1) mortality rates, 2) permanent ces- 
sation, both before first treatment and 3) heroin users 
that started their first treatment before the observation 
period, we need to consider possible limitations: 

1) Mortality rates were extrapolated by applying to the 
whole country figures from the North-East of Spain 



where heroin use injection was more frequent than in 
the rest of the country [22]. The extrapolation appears 
to be appropriate, since the period where the highest 
mortality rates are found for the two cohorts studied 
(1985 to 1999) coincides with the period when there 
were more HIV and drug injection related deaths in 
Spain [16]. However, if the extrapolation is not appropri- 
ate it would lead to over-estimation of the total inci- 
dence of heroin consumption for the whole country. 
Note that adding an additional 0.01 to the yearly mortal- 
ity rates, i.e. to account for the risk of dying when out of 
treatment being greater, would lead to even greater over- 
estimations of the incidence. 

2) Using lasting cessation rates from long-term cohort 
studies would overestimate incidence as they include per- 
sons with long cessation periods who finally may relapse. 
On the other hand, the fact that experimental users were 
not included in studies estimating cessation rates would 
produce underestimates. Nevertheless, these experimen- 
ters are only of anecdotic value for policy interventions. 

3) Although the exact dates and figures we have taken 
for first treatment probabilities prior to the observation 
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period may not apply to the whole country, the values 
assumed seem plausible as we obtained an increasing se- 
quence of probabilities from 1982 up to the first one 
estimated based on observed data in 1991, which hap- 
pens to have a similar slope to that estimated in the 
observed period. 

Besides model building, it is important to consider 
other limitations related to treatment data, both to its 
overall availability and to its accuracy. Treatment regis- 
ter data covers public and publicly funded centres, 
missing people using private treatment centres. This 
entails a small proportion especially once public substi- 
tution treatment centres were widely implemented all 
over the country following legislation in 1990 [28]. In 
relation to treatment variables used, we acknowledge 
the possibility of error in the reported year of heroin 
use onset, in which we cannot discern any systematic 
trend except perhaps a certain propensity to round to 
years ending in 0 or 5. 

Incidence trends by route of administration do not ne- 
cessarily reflect the route used at the time of onset, as 
the variable was collected referring to the 30 days prior 



to first treatment. However, in a previous study involving 
heroin users, both in and out of treatment, and a mean 
length of use of 10 years, more than 50% did not change 
their initial route of heroin administration [29]. Thus 
the study of incidence trends by route of administration 
in the period immediately previous to first treatment can 
provide an idea of the different patterns of heroin ad- 
ministration during the heroin epidemic in Spain [30]. 
The higher probability of entering treatment among 
individuals declaring injection in the previous month 
may possibly be related to a change to a more harmful 
route of administration. 

Conclusions 

With a simple conceptual model of heroin users' trajec- 
tories related to treatment demand, it has been possible 
to obtain approximations of heroin use incidence trends. 
Moreover, different assumptions made do not systemat- 
ically skew the conclusions. However, enhancing accur- 
acy of drug users' trajectories and an updating of new 
treatment admissions will further contribute to better in- 
cidence estimates. 
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